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1.0  INTRODUCTION  AND  PROBLEM  STATEMENT 

Analytical  techniques  are  needed  by  the  Air  Force  for 
predicting  stress  intensity  factors  of  cracked  metallic 
structures  subjected  to  rapid  thermal  pulses.  A  schematic 
of  such  a  problem  application  is  illustrated  in  Figure  1-1. 
General  purpose  numerical  techniques  such  as  finite  element 
methods  are  currently  available  for  the  solution  of  such 
problems.  However,  they  require  time-consuming  finite 
element  modeling  of  the  structural  configuration,  including 
the  crack,  and  detailed  thermal/stress  analysis  must  be 
performed  for  each  thermal  transient  to  be  addressed. 

The  modeling  is  further  complicated  by  the  extremely  steep 
stress  gradients  which  exist  in  the  vicinity  of  the  crack 
tip.  The  stress  singularities  at  the  crack  tip  dictate  the 
use  of  an  extremely  fine  finite  element  mesh  in  this  region, 
or  a  special  element  which  has  the  appropriate  crack  tip 
singularity  built  into  the  element  interpolation  functions. 
Furthermore,  multiple  finite  element  models  and  analyses  are 
required  for  a  single  structure  if  we  desire  the  stress 
intensity  factors  as  a  function  of  crack  size,  which  is 
usually  the  case. 

Thus,  the  objective  of  this  study  is  to  provide  a  convenient 
and  accurate  way  to  calculate  the  stress  intensity  factors 
caused  by  rapid  thermal  pulses,  for  use  in  conjunction  with 
existing  fracture  mechanics  software  used  by  the  Air  Force 
to  predict  crack  growth  and  fracture  in  flaw  tolerant  design 
applications. 

As  a  result  of  the  Phase  I  effort  reported  here,  a 
preliminary  version  of  an  IBM-PC  based  computer  program, 
AF-CRACK,  has  been  developed  based  on  the  concepts  of 


Green's  function  and  influence  function.  AF-CRACK 
calculates  stresses  and  stress  intensity  factors  by 
integrating  the  closed-form  Green's  function  solutions  to 
the  general  thermoelasticity  problem  of  a  point  heat  source 
in  a  flat  plate  containing  a  crack.  For  a  given  plate/crack 
geometry,  the  stress  distribution  and  stress  intensity 
factors  due  to  any  number  of  randomly  located  heat  sources 
(or  sinks)  can  be  calculated  by  AF-CRACK  in  a  few  minutes. 
Preliminary  verification  of  the  program  shows  that  results 
generated  by  AF-CRACK  are  very  close  to  those  obtained  by 
finite  element  methods. 

AF-CRACK  combines  and  links  several  independent  program 
modules  using  a  unique  "software  bus"  concept,  and  is  a  menu 
driven,  user-friendly  software  package  with  extensive 
graphics  capabilities.  Both  input  and  graphics  in  AF-CRACK 
use  the  popular  spread-sheet  program,  LOTUS-123  [1],  linked 
directly  into  the  AF-CRACK  program  through  the  software  bus. 
In  this  manner,  data  input,  manipulation  and  review  are 
greatly  facilitated  using  LOTUS-123 's  extensive  data 
management  and  graphics  capabilities. 

The  stress  intensity  factor  results  generated  by  AF-CRACK 
are  useful  as  input  to  any  of  a  number  of  fatigue  crack 
growth  and  fracture  prediction  programs  used  by  the  Air 
Force.  An  example  of  such  a  program,  Structural  Integrity 
Associates'  pc-CRACK  computer  program  has  been  included  on 
the  AF-CRACK  software  bus.  Although  a  direct  data  link 
between  the  two  programs  could  not  be  developed  under  the 
cost  and  schedule  limitations  of  the  current  Phase  I  effort, 
stress  intensity  results  generated  with  AF-CRACK  can  be 
directly  input  to  pc-CRACK  using  the  Kj  input  option  in  the 
pc-CRACK  LEFM  module.  If  a  Phase  II  effort  is  funded,  a 
direct  link  can  be  developed  between  AF-CRACK  and  any 
companion  fracture  mechanics  software  selected  by  the  Air 
Force . 


2.0  THEORETICAL  SOLUTION 


2 . 1  Assumptions 

The  following  assumptions  have  to  be  made  to  simplify  the 
problem: 

(a)  The  problem  is  assumed  to  be  plane  stress.  In  other 
words,  the  thickness  D  in  Figure  1-1  is  assumed  to  be  small 
enough  so  that  the  temperature  distribution  in  the  thickness 
direction  can  be  treated  as  uniform. 

(b)  The  metallic  materials  are  assumed  to  be  isotropic, 
homogeneous,  and  linear  elastic. 

(c)  The  rate  of  heat  application  is  slow  enough  that  the 
coupling  terms  and  inertia  terms  in  the  general 
thermoelasticity  equations  can  be  neglected,  i.e.  it  is 
assumed  that  quasi-static  thermoelasticity  applies.  It  is 
anticipated,  however,  that  Green’s  functions  could  be  found 
which  would  permit  the  incorporation  of  inertia  terms  into 
advanced  versions  of  the  software,  if  it  is  necessary  for 
envisioned  applications. 

(d)  The  crack  surfaces  are  assumed  to  be  stress  free  and 
heat  conductant.  Although  the  solution  techniques  developed 
in  this  report  can  also  be  extended  to  cracks  with  insulated 
surfaces,  for  purposes  of  this  Phase  I  feasibility  study, 
the  crack  surfaces  were  assumed  to  be  fully  heat  conductant. 
In  reality,  actual  crack  surfaces  are  expected  to  be 
somewhere  between  100  percent  heat  conductant  and  insulated. 

(e)  We  assume  that  the  heat  convection  coefficients  remain 
constant.  Again,  this  assumption  was  included  for  purposes 
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of  the  Phase  I  study.  Variable  heat  convection  coefficients 
could  be  incorporated  in  the  program  through  a  more 
complicated  but  similar  derivation. 

2.2  Governing  Stress  Solutions 


For  a  plane  stress  thermoelasticity  problem  as  illustrated 
in  Figure  1-1,  the  governing  equations  are 


fv2T  =  T,t  +  frj2T 


a  .  ■  .  =  0 


a..  =  2G  [£..  +  -  a(f^)6..T] 

ij  L  lj  ' 1— u 1  i]  kk  '1-y'  lj  J 

£ . .  =  (U.  .+u.  . ) /2 
13  V  1/3  3/iJ/ 


v2U.  +  (t-^)U.  ..  =  2(~^)aT 

l  v 1— u '  3,31  ' 1— u ' 


where  T^fx^x  ,t)  is  temperature  distribution,  t  and  x^  are 
time  and  Cartesian  coordinates  respectively,  f=(K/pc)  ,  K  is 
heat  conduction  coefficients,  p  is  mass  density,  c  is  heat 
2 

capacity,  77  =(2H)/K/D,  H  is  heat  convection  from  the  plate 
surfaces  to  the  environment  (see  Figure  1-1) ,  D  is  the  plate 
thickness,  are  displacements,  a  is  coefficient  of  thermal 
expansion,  and  G  and  u  are  shear  modulus  and  Poisson's  ratio 
respectively. 


Nominal  stress  distributions  generated  with  these  equations 
are  used  to  develop  stress  intensity  factors  via  the 
influence  function  approach  described  below. 

2.3  Stress  Intensity  Factor  Influence  Functions 

For  any  plate/crack  geometry,  [2]  cracktip  stress  intensity 
factors  can  be  determined  by  integrating  the  product  of  the 
stresses  at  the  crack  location  in  the  uncracked  structure 
and  an  influence  function  (or  weight  function) .  That  is, 
for  a  crack  as  shown  in  Figure  2-1,  the  stress  intensity 
factor  can  be  calculated  by 


K 


I 


K 


II 


CTyy(x)m1(x)dx 


axy(x)m2(x)dx 


(6) 

(7) 


where  m1(x)  and  m2(x)  are  influence  functions  and  cr^  are 
the  normal  and  shear  stress  distributions  on  the  cracked 
surface  in  an  indentical  but  uncracked  plate  which  is  under 
the  same  temperature  distribution  as  the  cracked  plate. 

Appendix  A  includes  influence  functions  for  a  single  edge 
cracked  plate  and  a  center  cracked  plate,  presented  by 
Bueckner  [2]  and  Tada  [3]  respectively.  Many  such  influence 
functions  for  other  crack  models  of  interest  are  available 
in  the  literature,  and  can  be  added  as  the  need  arises. 

2.4  Principle  of  Superposition 


In  general,  as  shown  in  Figure  2-2,  there  may  be  multiple 
heat  sources  (or  sinks)  within  the  plate  which  will  cause 


thermal  stresses  and  stress  intensity  factors  at  cracks  in 
the  structure.  Since  the  problem  is  linear,  it  is  easily 
demonstrated  that  total  stress  intensity  factor  can  be 
calculated  as  the  sum  of  the  stress  intensity  factors  caused 
by  each  individual  heat  source,  acting  independently.  That 
is 


K (t)  =  2  [K(t)]i 

i 


(8) 


where  K(t)  is  the  total  stress  intensity  factor  and  [K(t)]^ 
is  the  stress  intensity  factor  due  to  the  ith  heat  source. 

2.5  Green's  Function  for  Time  Integration 

In  equation  (8),  the  stress  intensity  factor  [K(t)]^  caused 
by  each  individual  heat  source  can  be  solved  with  the 
concept  of  Green's  function  integration.  As  illustrated  in 
Figure  2-3,  the  stress  intensity  factor  due  to  a  heat  source 
Q ( t)  can  be  calculated  by 


[K(t)  ]i  =  J  Qi(r)  G(t-r)  dr 

where  G(t)  is  the  stress  intensity  factors  due 
function  <5  ( t)  heat  source  at  the  l  heat  source 


(9) 

to  a  Delta 
location. 


As  shown  in  Figure  2-3,  generally,  the  Green's  function  G(t) 
will  decay  and  approach  to  zero  after  a  decay  period  t^. 
Therefore,  the  integration  range  in  equation  (9)  can  be 
reduced  to  from  (t-t^)  to  t.  Such  a  reduction  in  the 
integration  range  greatly  increases  the  speed  of  the 


calculation  because,  instead  of  integrating  for  the  entire 
time  history,  it  is  only  nessary  to  integrate  backwards  from 
the  present  time  t  to  (t-t^)  as  follows 


[K(t) ] i  =  J  Qi(T)  G(t-r)  dr  (9a) 

t“td 


2.6  Single  Edge  Cracked  Plate  and  Center  Cracked  Plate 


For  Phase  I,  only  two  crack  models,  single  edge  cracked 
plates  (Figure  2-4)  and  center  cracked  plates  (Figure  2-5) 
are  considered.  We  also  assumed  that,  similar  to  the  heat 
convection  coefficient  H  on  the  plate  surfaces,  the  heat 
convection  coefficients  H1  and  H2  at  the  edges  of  the  plate 
remain  constant. 


Previous  sections  show  that  the  stress  intensity  factors  for 
both  crack  models  can  be  easily  obtained  by  equations 
(6  through  9)  if  the  stress  distribution,  resulting  from  a 
delta  function  heat  source  (or  sink)  Q5(t)  at  any  arbitrary 
location  (x',y')  in  the  uncracked  plate  (Figure  2-6),  can  be 
solved. 


To  solve  for  the  Green's  function,  the  following  boundary 
conditions  need  to  be  included: 


at  x  =0, 


KT  =  -H  T 
,  x  i 

a  —  a  =  0 
xx  xy 


(10) 

(11) 
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at  x=B , 


KT  =  H_T 
,  x  * 

cr  =a  =0 
xx  xy 


at  y  =  ±°°, 


T  =  o 


a  =  a  =  a  =0 
xx  yy  xy 


The  temperature  solution  to  the  uncracked  plate  problem, 
shown  in  Figure  2-6,  can  be  obtained  in  closed  form  by  the 
standard  method  of  separation  of  variables  [4].  The  result 


T(x,y,t)  = - ^ -  2  {Z  (x)  Z  (x«) 

2pc6>/jrft  n 


exp[-fa^t-gT72t-^~^-'-^-  ; 


where 


anK(Hi+H2  ) 
tan  a  B  =  „  - - 

n  kV  -  H.  H, 
n  12 


Z  fx)  =  (Ka  cosa  x  +  H.  sina  x)  Y„ 
n  n  n  l  n  n 


2  2  2 
2 ) 


_ n  2  ' _ 

<K2«J+hJ  )  [  B  (K  )  +KH2  ]  +KH1  (kV+h^  ) 


The  next  step  is  to  solve  for  the  stress  distribution  due  to 
the  temperature  field  equation  (16).  References  [5,6]  show 
that  one  particular  solution  to  the  thermal  stress  problem 
can  be  expressed  in  terms  of  a  stress  function  $  as  follows 
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where  (J)  is  a  stress  function,  and  <j>o  is  a  function  of  x,  y, 
and  t  such  that  4>  would  remain  finite  as  t  approaches 
infinity.  Substitution  of  equation  (16)  into  equation  (24) 
yields 


4>  = 


-?n2t 


(25) 


where 
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and  erfc(x)  is  a  complementary  error  function. 

In  general,  the  particular  solution  4>  shown  in  equations  (25 
through  27)  does  not  satisfy  the  stress  free  boundary, 
equations  (11  and  13)  ,  at  x=0  and  x=B.  A  complementary 
solution  tf,  which  is  an  analytic  function,  must  be  included 
to  make  the  two  edges  x=0  and  x=B  stress  free.  The  stresses 
are  then  calculated  by 


a  =  2G  (*-<t>) 


a  =  2G  (^-4>) 
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where  is  the  solution  of 
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Solution  to  equations  (32  through  34)  can  be  obtained  by  the 
methods  described  in  [7]  for  the  problem  of  an  infinite 
strip  of  plate  subjected  to  arbitrary  tractions  at  both 
edges  (as  illustrated  in  Figure  2-7) ,  and  is  summarized  in 
Appendix  B. 

Therefore,  at  the  cross  section  of  y=0,  the  stresses  caused 
by  a  Delta  function  heat  source  Q6(t)  at  (x^y1)  are 


"yy  *  °yy  +  2  (Pncos“nx  +  °nsinanx> 


°xy  ’  axy  +  1  <Rncosanx  +  snsinanx> 


where  a  and  a  are  stresses  due  to  the  complementary 

yy  *y 

stress  function  and  are  listed  in  Appendix  B, 
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The  Green's  functions  for  the  stress  intensity  factors  can, 
thereby,  be  calculated  by  substituting  equations  (35  and  36) 
into  equations  (6  and  7)  ,  and  the  total  stress  intensity 
factors  can  be  obtained  by  summing  up  equation  (9)  for  all 
heat  sources  or  sinks  in  the  plate. 

> 

Verification  of  the  above  methodology  is  provided  in  Section 
3.3,  and  its  application  to  two  sample  problems  is 
illustrated  in  Section  4.0. 


K(t)= 

J  0 


<T  ™ 


Figure  2-1.  Influence  Function  Concept  for  Calculating 
Stress  Intensity  Factor 
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Figure  2-2.  Concept  of  Superposition  for  Multiple  Point 
Heat  Sources 
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3.0  PHASE  I  DEMONSTRATION  SOFTWARE  (AF-CRACK) 

3.1  General  Program  Features 

The  theoretical  solution  technique  described  in  Section  2.0 
has  been  implemented  in  the  form  of  a  demonstration  software 
package  called  AF-CRACK.  Although  preliminary  in  nature, 
this  software  package  incorporates  an  interactive, 
menu-driven  format,  which  permits  its  use  without  extensive 
skills  on  the  computer  or  operating  system,  and  without 
constant  reference  to  a  user's  manual.  Also,  because  it  is 
modularized,  a  user  can  stop  at  any  point  during  a  problem 
application  and  come  back  later  to  finish  the  calculation  at 
any  time. 

The  program  runs  on  an  IBM  PC-XT  or  -AT  (or  compatible)  with 
minimum  of  640K  memory.  A  math  co-processor  card  is  also 
recommended  (although  not  mandatory)  to  enhance  program 
execution  speed.  All  of  the  AF-CRACK  input  and  output  are 
unit-transparent.  The  user  can  use  any  set  of  units  he  or 
she  wishes  as  long  as  all  the  data  are  consistently  in  the 
same  unit  system.  The  following  sections  describe  the 
structure  and  execution  of  the  AF-CRACK  program,  as  well  as 
verification  problems  run  to  compare  the  results  to  other 
solution  techniques. 

3.2  Software  Architecture 

The  basic  structure  of  the  AF-CRACK  computer  program  is 
illustrated  in  Figure  3-1.  The  program  consists  of  six 
independent  modules  and  one  data  base  shared  by  all  the 
modules.  A  "software  bus"  concept,  developed  independently 
by  Structural  Integrity  Associates  [8],  was  used  to  link  the 
program  modules  and  the  data  base.  The  software  bus  allows 
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users  to  enter,  execute  and  exit  each  program  module  with  a 
single  key  stroke,  as  well  as  to  execute  other,  independent 
software,  which  can  be  set  up  to  interface  with  the  same 
data  base.  (Some  customization  of  data  files  is  required, 
however,  in  order  for  AF-CRACK  to  automatically  transfer 
data  with  other  independent  software  packages  on  the  bus.) 

As  illustrated  in  Figure  3-1,  AF-CRACK  currently  contains 
six  modules:  Geometry  Input,  Material  Input,  Heat  Source 
Data  Input,  K  Calculation,  pc-CRACK,  and  Review.  The  last 
two  modules,  pc-CRACK  and  Review  actually  incorporate  two 
pre-existing  software  packages  which  are  useful  in  pre-  or 
post-processing  the  AF-CRACK  results.  The  Review  module 
loads  up  the  popular  LOTUS-123  spreadsheet  and  graphics 
software,  which  is  used  to  generate  tabular  input  to  some 
program  modules,  and  to  create  tables  and  graphical  displays 
of  the  results.  Translation  routines  are  included  which 
automatically  link  AF-CRACK  input  and  output  files  with 
LOTUS-123. 

PC-CRACK  refers  to  another  menu-driven  fracture  mechanics 
program  developed  independently  by  Structural  Integrity 
Associates  [9],  which  possesses  a  wide  variety  of  crack 
growth,  critical  flaw  size  and  elastic  plastic  fracture 
machanics  capabilities.  The  long  term  objective  is  to  be 
able  to  directly  feed  the  stress  intensity  factor  results 
from  AF-CRACK  into  PC-CRACK  (or  some  other  similar  program 
of  Air  Force's  choice)  for  further  analyses,  such  as  fatigue 
crack  propagation,  critical  crack  size,  and  corrosion  crack 
growth.  However,  at  the  present  stage,  the  data  base 
generated  by  AF-CRACK  is  not  compatible  with  the  input 
format  of  PC-CRACK.  Full  communication  between  AF-CRACK  and 
PC-CRACK  can  be  achieved  later  if  Phase  II  of  this  project 
is  pursued. 
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Copies  of  these  two  independent  software  packages,  along 
with  user's  manuals  are  included  with  this  report,  as  part 
of  the  Phase  I  deliverable.  A  detailed  description  of  the 
other  AF-CRACK  program  modules  and  the  data  base  is  given  in 
the  following  sections. 


3.2.1  Geometry  Input  Module 


As  shown  in  the  menu  for  this  module  illustrated  in  Figure 
3-2,  this  module  has  two  steps:  select  crack  model  and  input 
plate  dimensions.  In  the  first  step,  users  are  asked  to 
choose  between  a  single  edge  cracked  plate  and  a  center 
cracked  plate.  More  crack  models,  such  as  crack  emananting 
from  a  hole,  semi-elliptical  surface  crack  in  a  half  space, 
and  quater-circular  surface  crack  in  a  quarter  space,  can  be 
included  in  Phase  II  of  this  project.  Once  the  crack  model 
is  chosen,  the  user  can  execute  the  second  step  of  inputting 
geometric  dimensional  data.  For  the  two  crack  models 
currently  implemented,  only  plate  width  B  and  thickness  D 
are  needed.  The  format  for  input  of  the  data  is  self- 
explanatory. 


3.2.2  Material  Input 


As  illustrated  in  Figure  3-3,  this  module  consists  of  three 
major  steps:  input  heat  convection  constants,  H,  and  H2 , 
input  heat  conduction  constants,  K,  p,  and  c,  and  input 
elastic  constants,  E,  v ,  and  a.  Once  again,  the  input 
format  is  self  explanatory. 


3.2.3  Heat  Source  Data 


As  shown  in  the  menu  in  Figure  3-4,  there  are  three  oteps  in 
this  module.  The  first  step  asks  the  user  to  input  number 
of  heat  sources  (or  sinks) ,  number  of  eigen  values  to  be 
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included,  time  increment  At  to  be  used  in  the  computer 
simulation,  and  total  time  for  the  analysis.  The  maximum 
number  of  eigen  values  that  can  be  used  is  currently  set  to 
20.  A  higher  number  of  eigen  value  terms  will  result  in 
longer  computer  execution  in  the  Green's  function 
calculation.  It  is  recommended  that  users  perform  a  simple 
convergence  study  to  get  an  optimal  number  of  eigen  values. 
The  program  will  choose  an  optimal  At  for  numerical 
computation  based  on  the  material  constants  and  geometrical 
dimensions  of  the  problem.  The  optimal  time  increment  At  is 
set  at  1/30  of  the  decay  period  td  for  the  problem.  In  this 
preliminary  version  of  AF-CRACK,  the  total  time  has  to  be 
less  than  250 (At). 

The  next  step  in  this  module  is  to  input  heat  source 
locations  (x',y')  as  well  as  the  starting  and  ending  time 
for  each  heat  source  or  sink.  When  the  second  step  is 
chosen  in  this  module,  the  user  only  needs  to  input  numbers 
in  appropriate  cells  in  the  LOTUS-123  spread  sheet,  which  is 
loaded  automatically.  Upon  leaving  this  spread  sheet,  the 
input  data  are  saved  in  a  text  file  for  later  use. 

The  third  step  in  this  module  is  to  input  heat  source 
intensity  versus  time  for  each  heat  source.  Again,  a 
LOTUS-123  spread  sheet  will  appear  on  screen  and  the  user 
only  needs  to  input  numbers  in  the  cells.  Note  that,  in 
this  spread  sheet,  the  time  t  which  appears  in  the  first 
column  is  the  time  relative  to  the  starting  time  specified 
in  the  last  step.  In  other  words,  if  a  heat  source  Q(t) 
does  not  start  to  generate  heat  until  t  ,  the  user  can  input 
the  starting  time  as  t  in  the  last  step  and  input  Q(t-tQ) 
in  this  step.  In  this  spread  sheet,  the  time  increment  has 
been  set  automatically  to  be  the  optimal  At  calculated  in 
the  first  step  and  the  user  should  not  attemp  to  change  it. 
In  the  current  demonstration  version  of  AF-CRACK,  only  a 


total  of  250  steps  is  allowed  for  each  heat  source. 
Finally,  we  must  mention  that  the  heat  source  intensity 
should  be  positive  for  heat  sources  and  negative  for  heat 
sinks. 

3.2.4  Calculate  K 

Stress  distribution  and  stress  intensity  due  to  specified 
heat  sources  are  calculated  in  this  module.  As  illustrated 
in  Figure  3-5,  there  are  four  steps  in  this  module. 

The  first  step  is  to  input  crack  length  for  the  selected 
crack  model.  The  second  step  is  to  calculate  Green's 
functions  for  all  the  heat  sources.  The  third  is  to  modify 
heat  source  intensity  versus  time  sources  if  the  user 
desires.  (This  step  can  be  skipped  if  the  user  has  already 
defined  heat  source  intensity  versus  time  curves  and  does 
not  want  to  modify  them.)  The  last  step  is  to  calculate 
stress  intensity  factors  for  the  selected  crack  model  under 
the  specified  heat  sources.  Two  sets  of  stress  intensity 
factors  are  calculated  for  center  cracked  plates,  one  set 
for  each  crack  tip.  Only  one  set  of  stress  intensity 
factors  is  needed  for  single  edge  cracked  plates. 

3.2.5  Data  Base 

The  input  and  output  of  all  the  above  modules  are  saved  in 
several  data  files.  There  are  two  types  of  data  files  used 
in  the  program:  text  files  and  binary  files.  The  text  files 
are  tagged  with  " . PRN"  as  their  extension  and  the  bipary 
files  are  tagged  with  ".DAT”  as  their  extension.  For 
example,  "GEOM.PRN"  is  a  text  file  saved  after  the  geometry 
input  module  is  executed  and  "HEAT . DAT"  is  a  binary  file 
saved  after  the  heat  source  data  input  module  is  executed. 


3.2.6  Review 


This  module  is  used  to  display  the  input  and  output  data, 
either  in  tabular  or  graphics  format,  on  the  screen.  Hard 
copies  of  the  displays  can  also  be  obtained.  The 
commercially  available  LOTUS-123  spreadsheet/graphics 
software  package  is  used  throughout  this  module,  adding  a 
wide  range  of  input/output  options  to  AF-CRACK.  LOTUS-123 
compatible  files  are  automatically  generated  of  all  key 
AF-CRACK  data,  which  can  then  be  plotted,  printed,  or 
otherwise  manipulated,  using  LOTUS-123's  broad  range  of 
capabilities. 

To  facilitate  program  use  by  those  who  are  unfamilar  with 
LOTUS-123,  extensive  use  of  "Macro"  commands  has  been 
included  in  all  the  spread  sheets.  Thus  the  user  simply 

presses  one  or  two  keys  as  instructed  on  screen  to  see  a 

specific  plot  of  his  or  her  results.  This  review  module 
(Figure  3-6)  puts  all  the  data  in  five  different  spread 

sheets,  which  are  called  from  the  Review  module  menu.  The 
first  option  obtains  a  spreadsheet  of  stress  intensity 
versus  time  results;  the  second  spread  sheet  contains 

Green's  functions  for  each  heat  source.  The  third 
spreadsheet  is  used  to  plot  stress  distribution  in  the 
uncracked  plate  at  the  cross  section  of  crack,  and  the 
fourth  is  used  to  show  heat  source  intensity  versus  time 
curves.  The  last  spreadsheet  provides  a  display  of  general 
problem  information,  such  as  crack  model,  crack  length, 
material  constants,  etc.,  in  the  analysis. 

3.3  Program  Verification 

A  preliminary  verification  has  been  performed  for  AF-CRACK 
by  comparing  the  AF-CRACK  results  with  finite  element  and 
pc-CRACK  [9]  solutions. 
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The  verification  problem  analyzed  is  a  single  edge  cracked 
plate,  as  shown  in  Figure  2-4,  with  the  following  input: 

B  =  10  inches,  D  =  1  inch,  a  =  2  inches 

H  “  Hx  =  H2  =  0 

K  -  0.0002579  Btu/ (sec-F-in) 

p  =  0.09734  lb/in3,  c  =  2.178  Btu/(lb-F) 

6  —  6 

E  =  10x10  psi,  u  =  0.3,  a  =  12.73x10  in/in/F 

Q (t)  =  H (t)  Btu  at  x'  =  5  inches,  y'  =  0  inch. 

Eigen  values  included  =  20 

where  H(t)  is  a  Heaviside  step  function.  This  verification 
problem  represents  a  10  inch  wide  thin  aluminum  plate  with  a 

2- inch  long  edge  crack,  under  a  unit  step  heat  source 
applied  at  the  center  of  the  plate,  right  in  front  of  the 
crack  tip. 

Stress  and  stress  intensity  factor  Kj  for  the  verification 
problem,  calculated  by  AF-CRACK,  are  depicted  in  Figures  3-7 
and  3-8  respectively. 

For  the  same  problem,  a  finite  element  program,  FEM2D  [10] 
was  used  to  predict  the  steady  state  temperature  as  well  as 
thermal  stress  distributions  for  the  uncracked  plate.  The 
finite  element  mesh  used  is  shown  in  Figure  3-9,  and  the 
resulting  stress  distributions  are  also  plotted  in  Figure 

3- 7  along  with  the  AF-CRACK  results.  Due  to  the  symmetry 
conditions  at  y=0,  only  half  of  the  plate  was  modeled  in 


Figure  3-9.  A  total  of  88  eight-node  isoparametric  elements 
and  303  nodal  points  were  used  in  the  finite  element  model. 
Figure  3-7  shows  that  the  stresses  predicted  by  AF-CRACK  are 
very  close  to  that  calculated  by  finite  element.  The  stress 
distribution  obtained  from  the  finite  element  analysis  was 
then  input  into  pc-CRACK  [9]  to  calculate  the  stress 
intensity  factors.  The  steady-state  stress  intensity  factor 
predicted  by  pc-CRACK  is  shown  in  Figure  3-8  in  conjunction 
with  the  AF-CRACK  solutions.  Again,  very  good  agreement 
between  the  two  solutions  is  observed. 

Another  steady-state  stress  intensity  factor  solution  has 
also  been  obtained  by  actually  modeling  the  crack  in  the 
finite  element  analysis  and  calculating  the  energy  release 
rate  by  a  path-independent  line  integral,  J'.  Details  of 
and  its  related  path-independent  line  integrals  are 
discussed  in  Appendix  C  of  this  report.  The  finite  element 
mesh  used  in  the  calculation  is  the  same  as  that 

illustrated  in  Figure  3-9  except  that  part  of  the  symmetry 
line  in  the  model  was  set  to  stress  free  to  reflect  the 
existence  of  the  crack.  The  energy  release  rate  J  =  is 
then  used  to  calculate  the  stress  intensity  factor  with  the 
relation  between  J  and  K: 


where  E  is  the  Young's  modulus  of  the  material.  Stress 
intensity  factor  predicted  by  the  line  integral  J|  is  also 
plotted  in  Figure  3-8  which  illustrates  that  the  steady- 
state  stress  intensity  factors  predicted  by  AF-CRACK  and  the 
other  two  independent  methods  are  within  1%  of  each  other. 
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More  verification  problems  are  needed  to  fully  test  the 
program  in  the  next  phase  of  the  project.  Numerically,  it 
is  necessary  to  run  more  finite  element  analyses,  both 
steady  r tate  and  transient,  for  more  crack  geometries  and 
sizes.  It  would  also  be  desirable  to  verify  the  program 
experimentally,  by  conducting  some  tests  to  measure  the 
stress  and  stress  intensity  factors  in  plates  or  other 
structures  under  rapid  thermal  pulses.  Verification  by 
testing  will  also  be  proposed  as  part  of  the  Phase  II 
effort. 

Nonetheless,  based  on  the  preliminary  verification  performed 
here,  we  conclude  that  the  AF-CRACK  methodology  is  capable 
of  predicting  stress  intensity  factors  for  structures  under 
rapid  thermal  pulses,  accurately,  and  with  only  a  small 
fraction  of  the  computer  time  and  manpower  required  by  the 
conventional  finite  element  analyses.  Full  development  of 
the  AF-CRACK  program  will  thus  greatly  enhance  the 
capability  of  the  Air  Force  for  predicting  stress  intensity 
factors  and  performing  subsequent  fracture  mechanics 
analyses . 
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Figure  3-1.  AF-CRACK  Software  Architecture 
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Figure  3-2.  Geometry  Input  Module  Menu 
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Figure  3-3. 
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USER  OPTIONS: 


+  1-  INPUT  CRACK  LENGTH 

2-  CALCULATE  GREEN"S  FUNCTIONS 
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4-  CALCULATE  X 

5-  TERMINATE 


+  Current  Line 
t  Moves  Up 
1  Moves  Doun 
4-1  Selects  Option 
Fx  Selects  Option  x 

x>10  Use  Shift 
*  Function  Key 
(Fll=  Shift+Fl) 


REV 1  EH  RESULTS  (REVH) 
MODULE  MENU 


USER  OPTIOHS: 

♦  1-  PLOT  K  VERSUS  TIME 

2-  PLOT  GREEN'S  FUNCTIONS 
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6-  TERMINATION 


*  Current  Line 
t  Moves  Up 
4  Moves  Down 
4-1  Selects  Option 
Fx  Selects  Option  x 
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4.0  SAMPLE  PROBLEMS 


Two  sample  problems  are  presented  in  this  section  to 
illustrate  the  use  of  the  AF-CRACK  program.  Except  for 
crack  model,  crack  size,  heat  transfer  coefficient  and  heat 
source  specification,  all  other  parameters  for  the  two 
sample  problems  are  the  same,  and  are  as  follows: 

B  =  10  inches,  D  =  0.2  inch 

K  =  0.0002579  Btu/ ( sec-F-in) 

P  =  0.09734  lb/in3,  c  =  2.178  Btu/(lb-F) 

E  =  10xl06  psi,  v  =  0.3,  a  =  12.73xl0~6  in/in/F 
Eigen  values  included  =  50 

The  above  material  properties  are  representative  of 
aluminium. 

4.1  Sample  Problem  1  -  Single  Edge  Cracked  Plate 

As  shown  in  Figure  4-1,  this  sample  problem  considers  a 
continuous  surface  crack,  2  inches  deep,  at  the  left  edge  of 
the  above  plate.  The  heat  convection  coefficients  used  in 
this  problem  are 


H  =  H1  =  H2  =  0.0001929  Btu/ ( sec-F-in2 ) 

and  the  four  heat  sources  are  at  (x,y)  coordinates  of 
(10,4),  (10,2),  (10,-2),  and  (10,-4)  with 
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Q4  =  H (t-150)  -  H ( t-1200 ) 


Btu 


Ql  =  Q2  =  Q3  = 


where  H(t)  is  a  Heaviside  step  function.  As  can  be  seen  in 
Figure  4-1,  this  simulates  a  line  heat  loading  at  the  right 
edge  of  the  plate,  which  steps  on  at  150  seconds,  and  off 
again  at  1200  seconds. 


Resulting  stress  distributions  for  the  uncracked  plate  are 
illustrated  in  Figure  4-2,  at  various  times  during  the 
loading.  The  stress  pattern  consists  of  compression  near 
the  heat  sources,  and  tension  near  the  edge  crack.  This 
shape  remains  constant  as  the  stress  builds  up  to  its 
maximum  value,  following  application  of  the  load,  and  then 
subsides  as  the  load  is  removed.  The  resulting  stress 
intensity  factors  versus  time  are  illustrated  in  Figure  4-3. 
From  this  figure,  we  can  see  that  the  mode  I  stress 
intensity  factor,  K^.,  begins  to  rise  at  the  time  of  load 
application  (150  seconds) ,  levels  out  at  a  maximum  value  at 
about  600  seconds,  and  then  decays  to  zero  again  after  the 
load  is  removed  at  1200  seconds.  As  expected  due  to  problem 
symmetry,  the  mode  II  stress  intensity  factor  is  shown  as 
zero  throughout  the  problem;  however,  this  was  not 
calculated  by  the  program  since  the  single  edge  crack  plate 
model  does  not  yet  incorporate  the  mode  II  influence 
functions.  Green's  functions  showing  the  stress  intensity 
factor  response  to  a  unit  spike  loading  at  each  of  the  heat 
sources  are  illustrated  in  Figure  4-4. 

4.2  Sample  Problem  2  -  Center  Cracked  Plate 

As  illustrated  in  Figure  4-5,  a  center  cracked  plate  with 
non-symmetric  heat  sources  was  analyzed  as  the  second  sample 
problem.  The  crack  length  in  this  problem  is  again  assumed 


to  be  2  inches.  Heat  convection  coefficients  used  are 


H1  =  H2  =  0,  H  =  0.0001929  Btu/ ( sec-F-in2 ) 
and  heat  sources  are  applied  at  (2.5,2)  and  (7. 5, -2)  with 

Ql  =  Q2  =  H (t-150)  -  H ( t-12 00 )  Btu 

As  can  be  seen  in  Figure  4-5,  this  corresponds  to  loading 
points  just  above  and  below  the  crack  plane,  on  alternating 
sides  of  the  crack,  and  subject  to  the  same  transient  load 
history  defined  above  in  sample  problem  1. 

Results  of  this  sample  problem  are  shown  in  Figures  4-6 
through  4-8.  Because  of  the  non-symmetry,  both  normal  and 
shear  stresses  develop  at  the  crack  plane,  as  illustrated  in 
Figures  4-6A  and  B.  The  stress  intensity  factor  shown  in 
Figure  4-7,  builds  up  and  decays  as  before,  but  in  this  case 
both  Kj  and  are  non-zero  and  are  calculated  directly  by 
the  program.  Figure  4-8  shows  the  Green's  functions  at  each 
crack  tip  (left  and  right)  which  result  from  the  unit  spike 
loading  at  heat  source  Q1> 


Qtt)(e  tu) 


/ 
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Figure  4-1.  Sample  Problem  I 
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Figure  4-2.  Stress  Distribution  at  y  =  0  for  Sample 

Problem  1  (a  =0) 
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Figure  4-6B.  Stress  Distributions  for  Sample  Problem  2  (a 
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STRESS  INTENSITY  FACTORS 


SAUPll  PROBLEM  2  (BOTH  CRACK  TIPS) 


5.0  CONCLUSIONS  AND  RECOMMENDATIONS 


The  Green's  function  concept  derived  in  this  report,  coupled 
with  the  well-known  use  of  influence  functions  for 
calculating  stress  intensity  factors,  is  shown  to  provide  a 
powerful  tool  for  fracture  mechanics  analysis  of  cracked 
metallic  structures  under  rapid  thermal  transients.  Due  to 
the  nature  of  the  Green's  function  decay  characteristics, 
stress  intensity  factors  for  selected  crack  models  can  be 
easily  calculated  in  a  few  minutes  on  an  IBM-PC  or 
compatible  personal  computer.  Such  a  short  turnaround  time 
makes  it  possible  to  conduct  thorough  parametric  studies  for 
cracked  structures  subject  to  this  type  of  loading,  at  a 
very  low  cost,  both  in  terms  of  computer  time  and  manpower, 
as  compared  to  other  means  of  problem  solution  known  to  the 
authors  at  this  time.  Also,  using  a  novel  software  bus 
concept,  the  resulting  computer  program,  AF-CRACK,  is  able 
to  interface  directly  with  the  popular  LOTUS-123 
spreadsheet/graphics  software  package,  which  greatly 
enhances  the  graphics  capability  and  ease  of  use  of  the 
program. 

Through  selected  verification  problems  for  which  other 
solutions  are  available,  the  AF-CRACK  program  is  shown  to 
accurately  predict  stress  and  stress  intensity  factor,  in 
comparison  to  other  techniques.  This  preliminary 
verification,  however,  covered  only  a  small  portion  of  the 
general  program  capabilities,  and  further  verification  is 
recommended . 

Based  on  the  success  of  the  Phase  I  feasibility  study 
reported  here,  we  concluded  that  it  is  highly  feasible  to 
develop  a  fast,  accurate  and  easy  to  use  general  purpose 
computer  program  based  on  this  methodology  to  predict  stress 
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intensity  factors  for  a  wide  range  of  metallic  (or  other) 
structures  of  interest  to  the  Air  Force,  under  rapid  thermal 
pulses.  It  is  also  possible  to  connect  the  program  with 
other  fracture  mechanics  computer  programs,  to  use  the 
resulting  stress  intensity  data  in  crack  propagation  or 
critical  flaw  size  predictions. 

Thus  we  recommend  pursuing  a  Phase  II  effort  to  complete  the 
program  development  effort  begun  by  this  study.  The 
following  areas  have  been  identified  as  candidates  for 
further  study/development  in  the  subsequent  phases  of  the 
project : 

(1)  Provide  more  crack  models.  The  basic  methodology  can  be 
expanded  applicable  to  a  wide  range  of  models,  such  as 
cracks  emanating  from  a  hole  in  a  plate,  cracks  in 
cylindrical  structures,  and  finite  aspect  ratio  cracks  (e.g. 
thumbnail  surface  cracks) . 

(2)  Extend  the  program  to  mixed  mode  crack  problems.  In  the 
preliminary  version  of  AF-CRACK,  only  one  of  the  two  crack 
models  has  both  mode  I  and  mode  II  solutions.  To  be  able  to 
handle  realistic  problems,  we  must  have  the  mixed  mode 
solutions  for  all  the  crack  models. 

(3)  Extend  the  Green's  function  solutions  to  problems  in 
which  the  crack  surfaces  are  insulated,  because  actual  crack 
surfaces  are  expected  to  behave  somewhere  between  fully  heat 
conductant  (as  assumed  here)  and  insulated.  Conceptually, 
this  type  of  Green's  function  can  be  derived  by 
superimposing  a  continuous  heat  dipole  along  the  crack 
surface  onto  the  existing  Green's  function 

(4)  Study  the  limitations  of  the  quasi-static 
thermoelasticity  used  in  this  study  relative  to  planned 


applications  (heating  rates)  of  the  software.  The  two 
assumptions  should  be  considered  independently  to  determine 
separately  at  what  heat  rates  the  lack  of  inertia  terms  and 
the  decoupling  of  the  equations  break  down. 

(5)  If  possible,  and  if  the  study  of  step  (4)  deems  it  to  be 
desirable,  include  inertia  effects  into  the  existing, 
decoupled  thermoelasticity  equations  for  the  Green's 
functions. 

(6)  Investigate  the  limitations  of  the  linear  elastic 
fracture  mechanics  theory  used  by  the  program  relative  to 
planned  applications. 

(7)  Perform  additional  verification  of  the  program  both 
through  additional  comparisons  with  existing  analytical 
solutions  and  through  experiments. 

(8)  Furnish  links  between  AF-CRACK  and  other  fracture 
mechanics  software  designated  by  the  Air  Force.  (eg. 
pc-CRACK,  CRACKS-84  &  -86,  CRACKGRO ,  ASDGRO ,  etc.) 

(9)  Conduct  parametric  studies,  using  the  software,  of  the 
effects  of  various  heat  sources  and  locations  on  stress 
intensity  factors  in  structural  configurations  of  interest 
to  the  Air  Force.  In  addition  to  providing  useful  technical 
results,  this  will  also  provide  feedback  to  the  program 
developers  on  user  friendliness,  speed  and  convenience  of 
the  program  in  a  typical  application. 
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APPENDIX  A 


INFLUENCE  FUNCTIONS 


A. 1  Single  Edge  Cracked  Plates 


For  a  single  edge  crack  in  an  infinite  strip  of  plate,  as 
illustrated  in  Figure  2-4,  the  influence  function  (or  weight 
function)  for  mode  I  cracking  has  been  provided  by  Bueckner 
[2] 
as 


m^x) 


[1 


,a-x, 
Pi<— > 


,a-x. 2 , 

p2(— )  i 


(Al) 


where 


p  =  0.6147  +  17.1844  R2  +  8.8722  R6 


(A2) 


P3  =  0.2502  + 


3.2889  R2  +  70.0444  R6 


(A3) 


R  =  a/B  (A4 ) 

and  x,  a,  and  B  are  defined  in  Figure  2-4. 

Bueckner  [2]  has  stated  that  the  above  equations  for  the 
weight  function  m^x)  are  very  accurate  for  cracks  up  to 
half  of  the  plate  width,  B.  The  authors  cannot  find  the 
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corresponding  influence  function  m2(x)  for  the  mode  II 
cracks  after  a  preliminary  literature  survey.  However,  we 
felt  that  the  solution  for  m2(x)  might  already  exist  or  can 
be  obtained  numerically  by  finite  element  analyses  in  the 
following  phase  of  this  project. 

A. 2  Center  Cracked  Plates 

Influence  functions  for  a  center  crack  in  an  infinite  strip 
of  plate,  as  illustrated  in  Figure  2-5,  have  been  derived  by 
Tada,  et  al .  [3]  as  follows 


m1(x) 


where 


Vx>  -  ^=<FiC|' -  T>1)(F2[|,(|  -  !>]> 


(A5) 


F1(f,rj)  =  1  +  0.297  J  (1  -  r,2 )  (1  -  cos|f ) 


( A6 ) 


1  ±  (sin^n/sin |f) 

■  --  = - -  (A7) 

,  1  -  (cos|f /cos~n)  2 

and  x,  a,  and  B  are  defined  in  Figure  2-5.  In  equation  (A7) 
the  plus  sign  is  for  the  right  crack  tip  and  the  minus  sign 
is  for  the  left  crack  tip. 

Tada  [3]  shows  that  equations  (A6)  and  (A7)  are  accurate 
within  1%  for  all  crack  siEes. 


F2(f  ,tj) 


- 


tanff 


A2 


APPENDIX  B  -  SOLUTION  FOR  COMPLEMENTARY  STRESS  FUNCTION  P 

This  Appendix  discusses  the  solution  procedure  for  the 

A  A 

stress  function  P  and  its  resulting  stress  a  and  a  in 

*  yy  *y 

equations  (28-36) . 

It  has  been  discussed  by  Timoshenko  [7]  that  the  solution 
for  the  problem  illustrated  in  Figure  2-7  can  be  written  as 


P  =  ^  J°°  (C1coshpx+C2sinh|3x+C3x  cosh/3x+C4x  sinh/3x)  cos/3y 


d/3 


+  ^  J0”*  (D1cosh/3x+D2sinh/3x+D3x  cosh/3x+D4x  sinh/3x)  sin/3y  d/3 


(Bl) 

where  x  =  x-(B/2).  In  equation  (Bl) ,  the  eight  constants 
Cl'  C2'  "  ‘  '  D3'  anc*  D4  are  determined  by  the  boundary 
conditions  at  x=0  and  x=B.  As  explained  in  equations  (32 
through  34) ,  the  boundary  conditions  for  P  at  the  two  edges, 
x=0  and  x=B,  are  the  non-zero  stresses  due  to  stress 
function  <|>  of  equation  (25).  By  substituting  equations  (25 
through  27)  into  equations  (33  and  34) ,  the  boundary 
conditions  for  P  can  be  written  as  follows 


at  x=0  ( x=  -B/2 ) , 


B2 


and  at  x=B  (x=B/2) , 


*.yy  "  VF> 
* , xy  -  F4<y> 


where 
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a ( 1+ u ) Q 
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exp(-w^) 
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V  irft 
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V- ir  a  Ft 
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2 

exp(-u  )  cj  exp(-w 
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V  irft  anft 
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exp(-u^) 

exp(a  (y-y'))  [erfc(u .)  +  -  — -  ] 

V  ttF  t  a 
b  n 


exp ( -w? ) 

-  exp (-a  (y-y'))  [erfc(u?)  +  —  - -  ])  (E 

n  V  irft  a 

*  n 


a(l+u)Q 

InrD 


F,(y) 


2  Z  ( B)  Z  (x ' )  { 


tv;  tv> 


exp(a  (y-y  • ) )  [a  erfc(u  )  +  2 


2  2 
exp(-(J1)  u1exp(-u1) 

V  irft  V~ ir  anft 


+  exp(-an(y-y ' )  )  [anerfc(u2)  +  2 


exp(-u2) 


u2exp(-(J2) 
V  rr  a  ft 


F4(y) 


2  Z*  (B)  Z  (X*  )  ( 

4pcD  nv  '  nv  7  1 


exp(-w  ) 

exp(a  (y-y'))  [-erfc(u  ) - - -  ] 

n  A  V  irft  a 

k  n 


exp(-w  ) 

exp  (-a  (y-y'))  [erfc(u?)  +  -  ]} 

n  2  an 


After  solving  the  eight  constants,  C^,  C ,  D3 ,  and  D4 , 
the  stresses  at  y=0,  can  be  expressed  in  a  double 
integration  as  follows 


pOO  pOO 

=  -  —  {[(32(R1  coshPx  +  R  sinh  Px 

"  J  -00  J  0  1  2 


+  R^x  coshPx 


+  R4x  sinhPx)  F^y)]  +  [P  ( — R ^  coshPx  +  R2  sinh  px 
+  R^x  coshPx-  R4x  sinhPx)  (y)  ]  )  cosPy  dp  dy 

p<»  pOO 

-  —  { [P2 (V  coshPx  +  v  sinh  Px  +  V  x  coshPx 


x  +  V3x  coshPx 


+  V  v  sinhpx) )  [F, (y)+F, (y) )  }  sinPy  dpdy 


(BIO) 


—  J  J  {  [  (V^p  smhPx  +  V2P  cosh  px  +  V3(cosh/3x 

—oo  Q 

+  Px  sinh/3x)  +  V4(sihPx  +  Px  cosh0x) ] 
[F3(y)+F4(y) ]  }  p  sinpy  dp  dy 


{  [ (R^/3  sinh/3x  +  R2P  cosh  Px 


-00  ^  o 


+  R3(cosh/3x  +  Px  sinh/3x)  +  R4(sih Px  +  Px  coshpx)  ]F  1(y) 
+  [-R^P  sinh/3x  +  R2P  cosh  Px  +  R3  (coshpx  +  Px  sinh]3x) 
-R4(sih/3x  +  Px  coshpx)  ]  F3  (y)  }  p  sinJ3y  d/3  dy  (Bll) 


where 


R1  =  -|-[sinh(|^)  +  (|5)cosh(|^)  ]/ (pB  +  sinh/3B)  (B12) 

P 


R2=  — —  [cosh ( )  +  (|^)sinh(|5)  ]/(/3B  -  sinh/3B)  (B13) 


R3=  -[-jp  cosh (§2)  ]/(PB  -  sinh/33) 


(B14 ) 


[-—  sinh(|^)  ]/ (PB  +  sinh/3B) 


( B15 ) 


[t§~  sinh(^)  1/  (PB  +  sinh/3B) 


( B1 6 ) 


v2=  [jf-  cosh(^)  ]/(/3B  -  smhpB)  (B17) 

V3=  -  [-j-  sinh(p)  ]/(/3B  -  sinhpB)  (B13) 

V4=  [~W~  cosh(p)  ]/(PB  +  sinh^B)  (B19) 

In  numerical  calculation,  the  integration  in  equations  (BIO) 
and  (Bll)  can  be  evaluated  by  using  the  Gauss-Laguerre 
quadrature  and  Gauss-Hermite  quadrature  [11].  In  AF-CRACK, 
six  points  were  used  for  the  Gauss-Laguerre  quadrature,  and 
five  points  are  used  for  the  Gauss-Hermite  quadrature.  That 
is,  the  integration  can  be  carried  out  approximately  by 

p° o 

F ( z )  dz  =  I  W.  F ( z . )  exp ( z . )  (B20) 

J0 

where 


Z±  =  0.222846604179 
Z  2  =  1.188932101673 
Z3  =  2.992736326059 

z .  =  5.775143569105 
4 

z  r-  =  9.837467418383 
b 

z,  =  15.982873980602 


W  =  0.458964673950 
W2  =  0.417000830772 
W3  =  0.113373382074 

W.  =  0.103991974531x10' 
4 

Wc  =  0.261017202815X10 
b 

W,  =  0.898547906430x10 
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APPENDIX  C  -  PATH  INDEPENDENT  LINE  INTEGRALS 
FOR  STEADY-STATE,  TWO-DIMENSIONAL  THERMOELASTICITY 

ABSTRACT 

Three  path- independent  line  integrals  J^,  M’,  and  L^  are  derived 
for  the  steady-state,  two-dimensional  thermoelasticity.  These 
integrals  are  similar  to  the  JK,  M,  and  L3  presented  by  Knowles 
and  Sternberg  [1],  but  include  additional  terms  of  either  free 

* 

expansion  displacement  vector  u^  or  temperature  9  and  its  complex 
conjugate  Q  in  their  formulation.  These  new  line  integrals 
enable  us  to  avoid  the  undesirable  area  integration  [2]  when 
calculating  the  strain  energy  release  rate  for  crack  problems. 
Application  of  J^,  M'  ,  and  L^  is  demonstrated  througth  a  sample 
problem  of  a  constant  heat  flux  disturbed  by  a  finite  crack  in  an 
infinite  plate. 

1.0  INTRODUCTION 

Since  the  discovery  of  the  conservation  integrals  J^,  M,  and  L  by 
Knowles  and  Sternberg  [1],  several  similar  conservation  integrals 
have  been  introduced  for  the  thermoelasticity.  Gurtin  [3]  has 
proposed  a  line  integral  which  consists  of,  in  addition  to  the 
well-known  J^-integral,  three  more  terms  related  to  the 
temperature  field.  Unfortunately,  path-independence  of  the  line 
integral  proposed  by  Gurtin  relies  on  several  restrictions  which 
in  general  can  not  be  met,  e.g.,  one  of  the  restrictions  is  that 
the  temperature  distribution  has  to  be  symmetric  about  the  crack 
axis  and  equal  to  zero  on  the  crack  surfaces.  Aoki,  et  al.  [2] 

.A  A  A 

have  derived  another  set  of  path  independent  integrals  J^,  M,  L, 

A 

and  I  for  general  elastic-plastic  problems.  As  a  special  case 
for  the  two-dimensional  thermoelast icity ,  these  integrals  become 
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5k  -  Jr<Wnk-Tjuj,k)ds  -  JA0mjamjal8  dA 

"  =  Jr(WnrTkuk,j)xjds  -  J/jk6jkal<0,mV'8>dA 

£3  *  Jrej3k'<Wnj-Vm,j>xk  +  Tjuk|ds 

-  e._,a.  (2a.  6.  9  +  a  6  0  .  x,  )  dA 

JA  ]  3k  lv  km  jm  pq  pq  ,3  k' 


(Cl) 


( C2 ) 


( C3 ) 


where  j,  k,  m,  p,  and  q  =  1  or  2,  W  is  strain  energy  density,  a 
is  equivalent  coefficient  of  thermal  expansion  defined  in  (C14) , 
0  is  temperature  distribution,  T  is  a  intergration  contour 
enclosing  the  crack  tip,  A  is  the  area  bounded  by  T,  and  e. is 
an  alternate  tensor.  Although  the  above  three  integrals  are 
path- independent ,  addition  of  the  extra  area  integration  at  the 
end  of  (Cl),  (C2),  and  (C3)  have  ruined  one  of  the  merits  of  the 

path- independent  integrals.  For  isothermal  elastostatic 

problems,  the  three  line  integrals,  J^,  M,  and  L  introduced  by 
Knowles  and  Sternberg  [1]  allow  us  to  take  advantage  of  its 
path-independence  and  calculate  the  energy  release  rate  based  on 
solutions  far  away  from  the  singular  crack  tip.  However,  when 
calculating  the  energy  release  rate  by  (Cl) ,  (C2) ,  or  (C3) ,  a 

highly  accurate  stress  and  temperature  solution  near  the  crack 
tip  is  necessary  in  order  to  evaluate  the  area  integration 
accurately,  implying  higher  computer  cost  and  potential  arguments 
on  how  accurate  is  enough  for  the  solution  near  the  crack  tip. 


For  the  steady-state,  two-dimensional  thermoelasticity,  the 
present  study  eliminates  the  need  of  area  integrations  by 
introducing  three  path-independent  line  integrals  J',  M',  and  LI. 

K  3 


2 . 0  FORMULATION 


2.1  Governing  Equations 

Consider  a  plane  strain  or  plane  stress  deformation  field  under  a 
steady  state  temperature  distribution.  It  is  assumed  that  no 
body  force  or  distributed  heat  source  are  present.  Governing 
equations  for  the  steady-state,  two-dimensional  thermoelasticity 
problem  are 


“jk.k  -  0 

in  B 

(C4) 

a  ■ ,  =  6  . ,  Xfc  +  2 lie  ■ . 

jk  jk  mm  j  k 

“  5jk(3X+24)a0 

( C5 ) 

ejk  -  <uj,k  +uk,j>/2 

( C6 ) 

°jk"k  ■  T? 

on  S 

a 

(C7) 

on  S 

( C8 ) 

u 

0  .  .  =  0 
>  J  3 

in  B 

( C9 ) 

0 

CD 

II 

CD 

on  ST 

(CIO) 

0  _  =  Q 
,  n 

on  SQ 

(Cll) 

where  summation  convention  has  been  implicitly  used  for  repeated 
subscripts  with  j,  k  =  1,  2,  B  is  the  elastic  body  with  S  as  its 
boundary,  ,  Su,  ST,  and  are  subsets  of  the  boundary  S,  n  is 
the  normal  direction  of  S,  a  is  coefficient  of  thermal  expansion, 
and  X  and  p.  are  Lame's  elastic  constants  defined  by 


C3 


2(1+«1) 


in  which  E1=E  and  v  ^=v  for  plane  stress  problems,  and 
2 

E^E/Cl-u  ),  and  u1=u/(l-u)  for  plane  strain  problems  (E  and  v 
are  Young's  modulus  and  Poisson's  ratio  respectively). 


2.2  Conservation  Laws 

As  a  start,  let  us  assume  that  the  elastic  body  B  is  simply 
connected.  Extension  of  the  formulation  to  multiply  connected 
bodies  is  discussed  in  2.4.  The  area  integrations  in  (Cl) 
through  (C3)  can  be  converted  into  line  integrations  by 
introductin  a  new  variable  ft,  which  is  a  complex  conjugate 
function  of  0,  into  the  formulation.  0  and  f2  are  thus  related  by 
the  Cauchy-Riemann  equations 


0  =  Q  ,  and  0  _  =  -ft  , 

fi-  ,  Z  ,  Z  ,  1 


(C13) 


Let  us  define  another  new  variable  u^ ,  the  free  expansion 
displacement  ,  as 


f* 

*  * 

U1  +  1U2  =  “l  J 


(0  +  iQ)dz 


(C14) 


where  z=x^+ix2  and  a^=a  for  plane  stress  problems  or  a1=(l+u)a 
for  plane  strain  problems.  With  (C13)  and  (C14),  the  area 
integration  in  (Cl)  can  be  rewritten  as 

a  n 

★  ★ 

a.  5.  a. 6  .  dA  =  a.  u.  ,  dA  =  a.  n  u.  ,  ds  (C15) 

J  jm  3m  1  ,k  J  jm  3,mk  J  jm  m  3  ,  k 


where 


Thus  (Cl)  becomes 


Jk  =  I  [Wnk  -  Vj.*1  ds 


where 


Uj  Uj  ~  Uj  ' 


or 


=  1 


{ Wn.  - 


T-u^  .  +  a  (T  0  +  T  Q) }  ds 


1  j  j,l  “l^lC 


and 


■J 


2  -  ,|  <Wn2  -  TjUj,2  -  “l<Tin  '  T2S>  I  d= 
r 


Similarly,  (C2)  and  (C3)  can  be  rewritten 


as 


(C17) 


(CIS) 


(C19) 


(C20) 


M'  =  Jr(Wnj-TkUk,j)Xjds 
or 


(C21) 


M' 


(  (Wnj-T^.lx.  +  a10(T1x1+T2>t2)-ain(T1x2-T2x1)  |d! 


and 


(C22) 


'3  *  Jrej3k|(Wnj-Tmu;,j)xk  +  Tjuk>ds 


(C23  ) 


L3  can  not  be  expressed  in  terms  of  0  and  n  like  (C19)  ,  (C20), 

and  (C22)  because  the  last  term  in  the  integrand  of  (C23)  would 
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result  terms  involving  /(0  +  if7)dz.  Thus,  for  the  steady-state, 
two-dimensional  thermoelasticity,  the  conservation  integrals 
defined  in  (Cl),  (C2) ,  and  (C3)  can  be  reduced  to  line  integrals 

with  the  addition  of  u*  or  0  and  Q  in  their  formulation.  We  see 
from  (C17 ) ,  (C21) ,  and  (C23)  that,  with  u^  being  the  parts  of  the 


displacement,  formulations  of  J£,  M' ,  and  for  the 
thermoelasticity  are  identical  to  formulation  of  J^,  M,  and  L3 
for  the  isothermal  elasticity.  This  observation  is  similar  to 
the  analogy  between  the  governing  equations  for  the 
thermoelasticity  and  isothermal  elasticity  discussed  in  [4]. 


One  way  to  calculate  J^,  M',  or  is  to  solve  for  D,  in  addition 
to  the  original  heat  transf er/thermal  stress  problem,  and  use 
( C19 ) ,  (C20)  or  (C21).  Since  Q  is  the  complex  conjugate  of  G,  it 
satisfies  the  Laplace  equation  and  can  be  solved  in  closed  form 
or  numerically  with  a  set  of  boundary  conditions  conjugate  to 
that  for  8.  Another  words,  Q  is  the  solution  of 


.  .  =  0 

JJ 

in  B 

(C24) 

n  =  0 

on  sT 

(C25) 

and 


n  =  -Q  on  S 

,  s 


(C26) 


where  n  and  s  are  normal  and  tangential  directions  respectively 
of  the  boundary  S.  Alternatively,  J£,  M',  and  can  be 

★ 


calculated  by  solving  for  the  free  expansion  displacement 
substituting  them  into  (C17),  (C21),  or  (C23). 


u  . 
J 


and 


> 


I 


C6 


2.3  Free  Expansion  Displacement  u. 


The  free  expansion  displacement  u ^  introduced  in  (C14)  is  not 
only  for  the  convenience  of  formulation  but  also  physically 
meaningful.  Substitution  of  (C16)  and  (C6)  into  (C5)  yields 


*  6jkx<uj,j+uk,k>  +  "<uj,k+uk,j> 


( C2  7 ) 


where  u_!  is  defined  by  (C18).  By  the  definition  of  u.!  ,  it  is 

★ 

easy  to  deduce  that,  when  uj=uj>  all  displacement,  stress,  and 

★ 

strain  components  vanish.  Therefore,  u^  is  the  stress  free 
displacement  for  the  elastic  body  under  the  same  temperature  0  as 

★ 

the  original  problem.  For  a  simply  connected  body,  u^  is  the 
displacement  solution  for  the  same  body  B  under  the  same 
temperature  distribution  0  but  with  homogeneous  mechanical 
boundary  conditions,  i.e.  T^=0  and  no  prescribed  displacements  on 
the  boundary.  It  is  worth  noting  that,  although  all  the  stress 

it 

components  corresponding  to  in  the  x^-x^  plane  are  zero,  the 

* 

stress  component  a  associated  with  u.  will,  in  general,  not  be 

J  J  J 

zero  for  plane  strain  problems. 


2.4  Multiply  Connected  Body 


Since  most  of  the  elastic  bodies  we  deal  with  in  fracture 
problems  are  multiply  connected,  it  is  important  to  extend  J.1, , 
M',  and  to  multiply  connected  bodies.  For  an  (m-rl)-ply 

★ 

connected  body,  Q  and  are  still  the  complex  conjugate  of  G  and 
the  stress  free  displacement  field  with  homogeneous  mechanical 
boundary  conditions,  respectively.  But,  when  solving  for  Q  and 

u*,  m  additional  cuts  have  to  be  introduced  to  make  the  body 
simply  connected.  Surfaces  of  these  m  cuts  are  stress  free  and 
can  move  freely  against  each  other,  i.e.  overlapping  or  sliding 
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of  the  two  adjacent  faces  on  each  cut  are  permitted  for  u* . 

Also,  fi  is  continuous  across  the  cuts  with  n  being  the  normal 
/  n 

direction  of  the  cuts.  The  additional  m  cuts  are  necessary 

because,  for  a  multiply  connected  body,  Q  and  u ^  may  be 
multi-valued  even  though  its  counter  part  0  is  still 

singled-valued  [4].  In  general,  Q  and  u*  are  not  continuous 
across  the  cuts,  in  a  multiply  connected  body,  therefore, 
formulations  for  J^,  M'  ,  and  must  be  modified  to  account  for 
the  discontinuities  when  the  integration  contour  T  intersects  any 
of  the  cuts.  For  instance,  for  contour  in  Figure  1,  an 

additional  term  of 


Tj[Uj,k]  ds 


(C23 ) 


VUj,lc1Xk  ds 


(C  29) 


and 


e  .  (T  [u*]  -T,  [u*  ]x.)  ds 

3mj '  mL  j  J  kL  k,mJ  j ' 


(C30) 


must  be  added  to  (C17),  (C21),  and  (C23)  respectively,  where  [f] 
is  the  discontinuity  of  f  across  C^.  Nevertheless,  it  is  always 
possible  to  chose  the  cuts  and  the  integration  contours  r  in  such 
a  way  that  they  do  not  intersect  each  other,  e.g.  contour  in 
Figure  1.  In  this  case,  expressions  in  (C17)  to  (C23)  are  still 
valid  for  multiply  connected  bodies  and  the  addition  of  (C28), 
(C29)  or  (C30)  to  the  conservation  integrals  can  be  avoided. 

Finally,  when  solving  for  n  and  u* ,  selection  of  the  m  cuts  to 
make  an  (m+l)-ply  connected  body  simply  connected  is  not  unique, 


neither  are  the  solutions  for  Q  and  u 


but  the  final  values  o! 


J^,  M',  and  do  not  depend  on  the  choice  of  the  cuts  or  the 

integration  contours. 

2.5  Energy  Release  Rates  and  Stress  Intensity  Factors 

Similar  to  the  physical  intepretat ion  of  J^,  M,  and  L3  in  the 
isothermal  elasticity  [5],  J£,  M',  and  defined  in  (C17)  to 

(C23)  can  be  deemed  as  the  energy  release  rates  associated  with 
translation,  rotation,  and  self-similar  expansion  of  the  crack, 
respectively.  Moreover,  since  the  order  of  stress  singularity  at 
the  crack  tip  for  the  two  dimensional  thermoelasticity  is  still 
-1/2  [6],  the  relationships  between  energy  release  rates  and 
stress  intensity  factors  for  the  isothermal  plane  elasticity 
[5,7,8]  also  apply  to  the  thermoelasticity  problems,  i.e. 

J[  =  (K^+KjI)/E1  ( C  3 1 ) 

and 

M'  =  x°J^  (C32) 

where  E1  is  defined  in  section  2.1  and  x°  is  the  x^coordmate  cf 
the  crack  tip  in  a  Cartesian  coordinate  system  of  which  x,-axis 
is  along  the  crack  surface. 

3 . 0  EXAMPLE  PROBLEM 

To  study  the  effectiveness  of  J'  M',  and  L'  in  numerical 

K  J 

calculation,  an  example  problem  of  a  finite  crack  of  length  2a  in 
an  infinite  plate  with  prescribed  temperature  gradient  vO  at 
infinity  is  chosen.  The  temperature  gradient  vO  at  infinity  is 
assumed  to  be  perpendicular  to  the  crack.  Stress  intensity 
factor  for  this  problem  has  been  found  [6]  to  be 


C9 


=  ( 2jiaaJ  ira  v0)/(1+k) 


(C33) 


where  p  is  shear  modulus,  a  is  coefficient  of  thermal  expansion, 
and  k  equals  to  (3-4u)  for  plane  strain  problems  or  (3-u)/(l+u) 
for  plane  stress  problems. 


A  finite  element  program  with  built-in  J  ^  routine  is  used  to 
solve  this  problem.  As  shown  in  Figure  C-2,  a  finite  size 
rectangular  plate  is  used  to  simulate  the  infinite  plate.  Due  to 
the  inherent  symmetry  of  the  problem,  only  one  half  of  the  plate 
is  needed  in  the  finite  element  analysis.  A  total  of  249  nodes 
and  70  eight-node  isoparametric  elements  is  used.  We  found  after 
a  few  numerical  experiments  that  numerical  solution  to  this 
problem  is  not  sensitive  to  the  size  of  the  finite  plate  once  the 
length  on  each  side  of  the  finite  plate  is  more  than  ten  times  of 
the  crack  length.  Also  shown  in  Figure  C-2  are  five  integration 
contours  to  calculate  J|.  A  standard  heat  transfer/thermal 
stress  analysis  is  first  performed  for  this  problem  followed  by 

it 

another  analysis  for  either  Q  or  u^  so  that  can  be  calculated 
at  the  five  seclected  contours  based  on  (C17)  or  (C20).  For  this 
doubly  connected  body,  a  cut  along  A-B  in  Figure  C-2  is 

★ 

introduced  when  solving  for  Q  and  u..  Since  there  is  no 

J 

intersection  between  the  cut  A-B  and  the  five  integration 
contours,  (C17)  and  (C20)  can  be  used  without  modification. 
Boundary  conditions  for  the  finite  element  analyses  are 
summarized  in  Table  C-l.  Once  J  ^  is  calculated,  the  stress 
intensity  factor  can  be  determined  by  (C31). 


Resulting  stress  intensity  factors  based  on  two  different 
equations,  (C17)  and  (C21)  are  listed  and  compared  with  the  exact 
solution  (C33)  in  Table  C-2.  It  is  seen  from  this  table  that 
stress  intensity  factors  predicted  by  J  ^  and  finite  element 
methods  are  essentially  path  independent  except  for  the  first 
contour.  Larger  errors  at  the  first  contour  is  expectable  since 


CIO 


the  finite  element  mesh  at  the  crack  tip  is  relatively  coarse, 
and  stress  distributions  near  the  crack  tip  are  thus  not  expected 
to  be  highly  accurate.  However,  with  only  249  nodal  points  and 
70  elements  in  the  finite  element  model,  stress  intensity  factors 
calculated  based  on  at  the  other  four  contours  are  within  5% 
or  3%,  depending  upon  whether  (C17)  or  (C21)  is  used,  of  the 
exact  solution  given  by  Sih  [6].  Judging  by  numerical  accuracy 
and  computational  cost,  formulations  of  J  '  ,  M',  and  L  '  in  terms 

K  J 

of  0  and  Q,  ( C19 )  ,  (C20)  and  (C22),  is  preferable  to  those  in 

★ 

terms  of  u ^  ,  (C17),  (C21)  and  (C23).  However,  formulations  of 

★ 

J',  M',  and  L'  in  terms  of  u.  still  provide  a  good  alternative 

K  -i  J 

when  n  is  difficult  to  obtain. 

4.0  CONCLUSION 


Conservation  integrals  for  thermoelasticity  have  been  found  to 
contain  extra  area  integration  terms  which  often  make  these 
path- independent  integrals  less  attractive  compared  to  other 
methods,  such  as  special  crack-tip  elements.  However,  for  plane 
strain  or  plane  stress  problems  under  steady  state  temperature 
distributions,  such  a  undesirable  area  integration  can  be 
eliminated  by  using  one  of  the  three  path- independent  line 


integrals,  J^,  M',  and  ,  introduced 

★ 

M  ,  or  ,  an  auxiliary  variable, 
addition  to  0  and  u^  of  the  original 

★ 

auxiliary  problem  of  u^  or  f?  can  be 
closed  form  without  much  difficulty. 


in  this  paper.  To  use  J,'., 

K, 

or  Q  ,  has  to  be  solved  in 
problem.  Solution  to  the 

obtained  numerically  or  in 
Physically,  n  is  a  complex 


conjugate  of  the  temperature  distribution  0  and  u*  are  thr  stress 
free  displacements.  It  has  been  shown  in  this  paper  through  an 
example  problem  that  the  energy  release  rates  or  stress  intensity 


factors  for  cracks  in  a  two  dimensional  solid  under  steady 
temperature  can  be  easily  calculated  with  J '  or  M' 
relatively  coarse  finite  element  model. 


state 


and 


a 


5. 
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Table  C-2 


Normalized  Stress  Intensity  Factor 


( l+(c )  K 


II 


2^aa'i  it  a  v0 


with  (2.11) 

with  (2.13) 

contour 

1 

1.133 

1.184 

contour 

2 

0.951 

0.974 

contour 

3 

0.955 

0.977 

contour 

4 

0.952 

0.975 

contour 

5 

0.952 

0.970 

Sih  [6] 
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1.000 


